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Abstract 

Molecular Dynamics is applied to Ferroelectric Perovskites in the framework of a first-principles 
derived effective Hamiltonian (Zhong, Vanderbilt, Rabe, Phys. Rev. Lett. 73 (1994), 1861). 
The degrees of freedom, that obey the Newton equations of motion, are the local modes and the 
displacement modes. The Nose-Hoover method is implemented, as well as the Parinello-Rahman 
scheme to perform fixed temperature and fixed stress tensor simulations. This allows to study the 
thermodynamics of ferroelectric perovskites and to reproduce successfully the Monte Carlo results 
on phase transitions, polarization and homogeneous strain evolution with temperature of BaTiOs, 
taken as an example. 
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I. INTRODUCTION 



The study of ferroelectric (FE) and related perovskite oxides has been the subject of many 
atomic-scale theoretical and numerical studies for at least fifteen years. In the ninetees, 
there has been an increasing growth of numerical work due to the appearance of significant 
computational facilities. Ab initio calculations have considerably developed, and have been 
used very early to derive parameters of phenomenological models used to understand and 
predict the thermodynamics of phase transitions in perovskite oxides, that can be quite 
complex. 

One of the most powerful tool, that captures most of the thermodynamics of FE mate- 
rials, is the so-called "Effective Hamiltonian" (EH), introduced by Zhong, Vanderbilt and 
Rabe in 1994. It can be viewed as a simplification of the potential energy surface of 
a ferroelectric in terms of relevant degrees of freedom such as the local modes, the inhomo- 
geneous strain tensor and the homogeneous strain tensor. The expression of this hamilto- 
nian is phenomenological and all the parameters are directly derived from first-principles 
DFT calculations, usually performed in the framework of the Local Density Approximation 
(LDA). Given this hamiltonian, one can perform Monte Carlo (MC) simulations to predict, 
at given temperature and pressure (or stress tensor), various physical quantities obtained 
as statistical averages (polarization, strain, dielectric constant, etc). Zhong, Vanderbilt and 
Rabe used it to reproduce successfully the complex sequence of phase transitions of barium 
titanatejl, 3|. 

The MC framework is a numerical technique that allows to sample the space of configu- 
rations according to an equilibrium probability distribution. An alternative technique, that 
normally yields similar results, is the Molecular Dynamics (MD). In MD, the degrees of 
freedom evolve with time according to the Newton equations of motion, leading, if the tra- 
jectory is long enough, to an equilibrium sampling of the phase space. Thus, time-averaged 
microscopic quantities on such equilibrium trajectories can be used to obtain macroscopic 
quantities, that should be the same as those obtained within MC simulations (ergodicity 
hypothesis). 

n 

MD has already been used with the Effective Hamiltonian. For instance, Burton et al[3\ 
applied it to the simulation of relaxor materials. Recently, Ponomareva et al applied it to 
the study of THz dielectric response of barium titanate J|. 



2 



In the present work, we apply the MD method to the Effective Hamiltonian. We imple- 
ment the Nose-Hoover and Parinello- Rahman schemes, that allow to perform MD with fixed 
temperature and fixed stress tensor. We first recall in the theoretical background the basis 
of the Effective Hamiltonian introduced by Zhong, Vanderbilt and Rabe, and describe the 
Nose-Hoover and Parinello-Rahman algorithms. Then the method is successfully applied to 
barium titanate to reproduce results very similar to those obtained by MC simulations. 



II. THEORETICAL BACKGROUND 



A. Degrees of Freedom 

The main idea of the effective hamiltonian is to decrease the number of degrees of freedom 
(a priori 15 per ABO3 unit cell in a perovskite) by using new coordinates representing 
collective motions of the atoms in each unit cell i. 

We assume that the considered perovskite oxides have an energy landscape which is well 
described in the framework of the effective hamiltonian proposed by Zhong et aim 1 994 [l], [3] - 
This hamiltonian deals with three kinds of degrees of freedom: 

(i) the local modes Ui, that roughly represent the polar displacements of the atoms of 
unit cell i, which induce an electric dipole in each unit cell. The dipole is obtained as Z*Ui, 
where Z* is the effective charge of the soft mode, calculated in the cubic structure. 

(ii) the inhomogeneous displacement modes Vi, that describe the local strain in cell i and 
refer to long-wavelength acoustic modes. From this displacement field, one constructs easily 
the inhomogeneous strain tensor field ^(z)^. 

(iii) the homogeneous strain tensor -qf, that does not depend on the unit cell and is a 
global deformation applied to the whole system. 

The strain existing at cell i is therefore T)i(i) = r}f + r]\{i). 

In many perovskite oxides, other atomic motions arising from the freezing of zone- 
boundary modes play a very important role. These acoustic modes at the R or M points of 
the First Brillouin Zone of the cubic structure, called antiferrodistortions (AFD), consist of 
rotations of the oxygen octahedra around a given axis. These modes are not considered here 
but could be dealt with exactly in the same way as local modes and displacement modes. In 
the following, we assume that the system studied consists of N unit cells. Thus the number 
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of internal degrees of freedom is 6N. 



B. Effective Hamiltonian 

The form proposed by Zhong et al in 1994|l,|2|] is : 



H e// ({u 4 },K«},{^})=^({^}) 
+E d v\{u l }) + E shoH {{u l }) 

+E das ({vt(i)} , {nF}) + , WS)} , {v?}), 

in which we can find a local mode self-energy E sel f (acting as a local potential for the 
local modes), a dipole-dipole electrostatic interaction energy E dpl , a short-range interaction 
energy E short (extended up to third-neighbors), an elastic energy E elas (depending only on 
the homogeneous and inhomogeneous strains) and a term coupling the strain to the local 
modes E mt (crucial to describe piezoelectric effects). The precise expression of all these 
terms is given in Ref. |2|. 

Interestingly, this hamiltonian only includes local anharmonicities in the onsite part and 
in the coupling between local modes and strain. 

C. Dynamics 

We assume that a mass can be affected to the local modes (mi m ) and displacement modes 
(m>dsp)- It is thus possible to formulate time-evolution equations for these quantities: 

where Im and dsp stand respectively for local mode and displacement mode. In the code 
we have implemented, the forces are calculated directly according to an analytical derivation 
of the effective hamiltonian (taken exactly as it is in Ref. Q). The masses mi m and m dsv are 
chosen according to Ref. |5. 
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III. MOLECULAR DYNAMICS 



The previous equations of motion are numerically integrated within the well-known Verlet 
algorithm^. A time step h is chosen, usually in the range 10~ 15 s. In this scheme, the local 
mode at time t+h is computed from the local modes at t and t — h, using the forces calculated 
at t: 

u^t + h) = 2u.it) -U i (t-h) + —fl m (t) (3) 
and idem for the displacement modes: 

Vi(t + h) = 2v.it) - Vi(t -h) + —ft p (t) (4) 

Of course, at each step, the different degrees of freedom are coupled and are coupled to 
the homogeneous strain: 

i~lra(±\ rl 



frit) = fi m ({uj(t)} , wax*)} , {fie®}) 



The dynamics is started from initial positions and from initial velocities, chosen randomly 
to reproduce the Maxwell distribution corresponding to an initial temperature, as usual in 
MD. 

For the first step, the positions and velocities are computed from a simple Taylor expan- 
sion: 

v^h) = v,(0) + /Ao) + -^/f p (0) 
at 2m dsp 

The instantaneous temperature is calculated from the velocities, based on equipartition 
theorem. 

2 



2 x ^Nk B T(t) = J2 \rn ln |^(*)} 

i ' 

+£^{^ (t) . 



The computation of the electrostatic dipole-dipole energy and electrostatic dipole-dipo 
forces is performed by the Ewald method, following the scheme proposed by Zhong et al^2\ 



e 
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A. Microcanonic molecular dynamics (NVE) 



In the microcanonic scheme, the homogeneous strain is fixed and does not vary all along 
the simulation: = rjf 1 ' . The temperature can not be chosen a priori and is found at 

the end of the simulation as the time average of the "instantaneous temperature". 



B. Nose-Hoover molecular dynamics (NVT) 

The Nose-Hoover algorithm is a method to fix the equilibrium temperature to the desired 
temperature T . We describe briefly this method in this section. A fictitious degree 

of freedom s with "mass" Q is added that represents a thermostat. The equations of motion 
are modified according to: 

m ds P ^p = fi ~ m dsp ((t) — 



d( _ k B gT(t) Tp 

dt Q [ T(t) ' 

where T(t) and T are respectively the instantaneous and targetted temperature, g is the 
number of degrees of freedom in the system, and ((t) = ^rr(t)- Q has not the dimension of a 
mass, but is the product of an energy by the square of a time. Its value must be appropriate 
so that the exchange of energy between the system and the thermostat are neither too slow 
nor too fast. Anyway, when chosen correctly, the precise value of Q does not influence the 
quality of the phase space sampling and has no effect on the macroscopic averages performed 
on an equilibrium trajectory. 

These equations are also numerically integrated within the Verlet algorithm: 

*(t + h) = 2u.it) - u^t -h) + —(fl m (t) - m lm C(t)^(t)) 

m lm dt 

in which the velocity at time step t is estimated according to[10]: 

dvii 3ui(t) - 4ui(t -h)+Ui(t-2h) , 

ir (t) = 2h (5) 
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C. Parinello- Rahman molecular dynamics 

The Parinello-Rahman algorithm[ll| is a very efficient method to perform molecular 
dynamics under fixed stress tensor. Such a method is crucial to study phase transitions in FE 
oxides since various phases currently appear (cubic, tetragonal, orthorhombic, rhomboedral, 
etc). In the Parinello-Rahman scheme, the supercell is allowed to vary in shape and volume 
along the simulation so that the macroscopic stress tensor (average of the instantaneous 
stress tensor) is equal to the desired one. We summarize the Parinello-Rahman method in 
the following, as we have implemented it. 

We denote by H(t) the 3x3 matrix formed by the components of the three vectors a, b and 
c that define the supercell of the simulation, expressed in an absolute orthonormal coordinate 
system, attached to the cubic 5-atom unit cell. It depends on t since the supercell is allowed 
to evolve in shape and volume along the simulation. H(t) is related to the homogeneous 
strain tensor r] H {t). 

The matrix G(t) is defined by G(t) = l H(t).H(t). The local mode at (i) has three 
components u' ia (a=x,y,z) in the time-dependent basis (a, b, c) and we denote by u[ the 
vector with components u' ia . We have Ui = i^(i).u-, where Ui is the vector formed by the 
components of the local mode at (i) in the absolute coordinate system. In the Parinello- 
Rahman framework, the equations of motion are formulated on the scaled variable u' ia : 

With the same notations for the displacement modes V{ = if (t).-^, 

The H(t) matrix evolves according to: 

r ^ r d?H . . . . 
W-^{t) = {a-a 1 ).u 

in which If is a "mass" associated to the dynamics of the unit cell vectors a, b and c, 
cjj - = dQ/dHij with Q the volume of the MD supercell, a is the instantaneous stress tensor 
and <7o the targetted stress tensor. The mass W has to be chosen correctly, in the same 
manner as Q. Its precise value, when chosen in the correct range, does not influence the 



thermal averages over equilibrium trajectories 
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Here again the equations of motion are integrated within the Verlet algorithm. The veloc- 
ities are calculated at each step according to Eq. [5l Using this algorithm, the instantaneous 
stress tensor components oscillate around the targetted value op so that their average 
over an equilibrium trajectory is equal to crp . 

The instantaneous stress tensor is calculated in this work by 

where Q is the supercell volume and the (i, j)-component of the homogeneous strain 
tensor. 



IV. EXAMPLE: BATIO3 

A code has been written in fortran 90, that includes microcanonic, Nose-Hoover, 
Parinello-Rahman and also both Nose-Hoover and Parinello-Rahman algorithms, that can 
be coupled to perform fixed temperature and fixed stress tensor MD^. This is precisely 
what is required to study the thermodynamics of ferroelectric materials and perform ther- 
mal averages. We present here the results obtained on the prototypical example of BaTiOs 
(BTO). The following computations have been performed with a time step of 2xl0~ 15 s on 
a supercell 12 x 12 x 12 with periodic boundary conditions. 

Barium titanate, the prototypical ferroelectric oxide, has been the first solid to which 



ler with Monte Carlo 
The local modes are 



the formalism of effective hamiltonian was applied in 1994, toget 
simulations. We use here the numerical parameters given in Ref. 
initialized to zero, as well as the displacement modes and the homogeneous strain. Each 
equilibrium trajectory is obtained independently from the others (i.e. the simulations are 
not initialized with the result of a previous computation). 

The simulations are performed in the Nose-Hoover and Parinello-Rahman scheme with 
an imposed diagonal stress tensor corresponding to a negative hydrostatic pressure P = -4.8 
GPa (i.e. positive stress tensor components). This negative pressure is the one used by 
Zhong et al in their simulations to correct the underestimation of the lattice constant by 
the LDA. We perform simulations of 10 5 steps. During the 50 000 first steps, the system 
is equilibrating, and the averages are computed on the 50 000 last steps. In fact, the 
thermodynamic equilibrium is in general quickly reached (within a few thousands steps). 
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We first give an example at T=210K, for which BTO is found orthorhombic. We examine 
the time evolution of the space-averaged homogeneous strain tensor components and of the 
space-averaged local modes (proportionnal to the instantaneous polarization). After « 4000 
steps, two components of the average local mode dissociate from the third and reach the 
value of ~ 0.024 ao, while the third one oscillates around zero. The 10 5 steps of this evolution 
are shown on Fig.CD Note that two states of macroscopic polarization (positive and negative) 
are of course possible for each component and degenerate. 



irii?i ripoLU |E30~| psGejartE-aaE-dg 



FIG. 1: Time evolution of the three components of the local modes (averaged over the whole 
simulation cell) at T=210K. 

It can be noted that the fluctuations of the vanishing component are, as expected, much 
larger than the two others. 

The time averaged local mode components (averaged over the whole simulation cell) and 
the time averaged homogeneous strain tensor components are drawn on Figs. [2] and [3] as 
a function of temperature. In these two figures, we have renamed the components so that 
the polarization keeps the same direction for the whole range of temperature of a given 
phase (since all our points are obtained independently from each other, there is absolutely 
no reason that, for example, all the tetragonal configurations have their polarization along 
z). As expected, we find the three phase transitions of BTO, well described as being first- 
order (both the polarization and the strain experience a discontinuity at the three phase 
transitions). 
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At low temperature, the three components of the macroscopic polarization are equal, as 
well as the three diagonal components of the homogeneous strain tensor. The non-diagonal 
components of the strain tensor are small and equal, but clearly non zero. This describes 
the rhombohedral (R) phase of barium titanate (space group R3m). 

At about 190 K, one of the components of the macroscopic polarization falls to zero. In 
the same way, one of the diagonal components of the homogeneous strain tensor dissociates 
from the two others and becomes smaller. This describes the orthorhombic phase of barium 
titanate (space group Amm2). We note that the non-diagonal components of the strain 
tensor are not all equal to zero: 774 (corresponding to 2e yz ) is not zero if rj2 (e w ) and 773 (e zz ) 
are the largest diagonal components. 

At about 222-229 K, another component of the macroscopic polarization falls to zero. 
The homogeneous strain tensor has one diagonal component larger than the two others, the 
non-diagonal ones being zero. This is the tetragonal phase of barium titanate (space group 
PAmm). 

Finally, at about 290-295 K (Curie temperature), all the components of the macroscopic 
polarization fall to zero, corresponding to the high-temperature cubic phase of barium ti- 
tanate (space group Pm3m). The diagonal components of the strain tensor are equal around 
1.012. Since the strain is referenced in this hamiltonian to the LDA ground state of BTO 
(a = 7.46 a.u.), this corresponds to a lattice constant of ~ 3.995 A. All the non-diagonal 
components of the strain tensor are zero. We note that our simulation does not evidence any 



thermal expansion, which is one of the well-known problems of the Effective Hamiltonian|l2|. 

The absolute values of polarization and homogeneous strains are very close to those of 
Zhong et a/fl Q, while the phase transition temperatures are slightly lower (~ 290-295, 
222-229 and 190 K in our case versus « 300, 230 and 200 in the case of Zhong et at). 

We now recalculate the temperature evolution of the polarization and strain of BaTiC>3, 
except that each temperature is now initialized with the result of a previous calculation. We 
perform these calculations by decreasing T from 360 to 160 K. 

Practically, the two last configurations (including local and displacement modes and also 
strain tensor) are stored at the end of each calculation, so that the Verlet algorithm can be 
restarted with a new temperature. The calculations are still performed in the Nose-Hoover 
and Parinello-Rahman framework, with an external pressure -4.8 GPa, corresponding to a 
positive stress tensor a xx = a yy = a zz =4.8 GPa (and a xy = a xz = a yz = 0). We perform 
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FIG. 2: Average local mode components as a function of temperature. R, O, T and C denote 
respectively the Rhombohedral, Orthorhombic, Tetragonal and Cubic phases of BaTiC^. Around 
phase transitions, on a small temperature range rs 10 K, either one phase or the other can be 
obtained since all the points are obtained independently from each other (The MD is in each case 
initialized in the same way). 
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FIG. 3: Average strain tensor components as a function of temperature (Voigt notation). R, O, 
T and C denote respectively the Rhombohedral, Orthorhombic, Tetragonal and Cubic phases of 
BaTiOs. Around phase transitions, on a small temperature range « 10 K, either one phase or the 
other can be obtained since all the points are obtained independently from each other (The MD is 
in each case initialized in the same way). 
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50000 steps for each temperature, with the same time step of 2x10 s. The macroscopic 
averages are computed on the 40000 last steps. 

As a consequence of this procedure, the temperature evolutions obtained are much more 
regular, especially close to phase transitions. Out of the regions close to phase transitions, 
the results are of course the same as in the previous case (within the numerical accuracy). 
The temperature evolution of polarization and strain are shown on Figs. [4] and [5] for the 
case of increasing temperature. The phase transitions are localized approximately at 190±5 
(R-O), 230±5 (O-T) and 295±5 K (T-C). 
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FIG. 4: Average local mode components as a function of temperature. R, O, T and C denote 
respectively the Rhombohedral, Orthorhombic, Tetragonal and Cubic phases of BaTiC>3. Results 
obtained by decreasing the temperature and initializing the MD from last two configurations of the 
previous temperature. 50000 steps are performed for each point, and the averages computed on the 
last 40000. 



V. CONCLUSION AND PERSPECTIVES 

The MD method with Nose-Hoover and Parinello- Rahman schemes has been applied suc- 
cessfully to barium titanate in the framework of the Effective HamiltonianQ. A code has 
been implemented in fortran 90. Thermal averages can be computed over the equilibrium 
trajectories obtained with these two techniques. This allows to reproduce the phase tran- 
sitions and the temperature evolution of the polarization and strain tensor components of 



barium titanate, as it is found from Monte Carlo simulations l[ 
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~\. The code implemented 
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FIG. 5: Average strain tensor components as a function of temperature (Voigt notation). R, O, 
T and C denote respectively the Rhombohedral, Orthorhombic, Tetragonal and Cubic phases of 
BaTiOs. Results obtained by decreasing the temperature and initializing the MD from last two 
configurations of the previous temperature. 50000 steps are performed for each point, and the 
averages computed on the last 40000. 

can thus a priori be used with other Effective Hamiltonian parameters that can be found in 
the litterature. 

The MD simulations can be used to get insight, not only into thermal averages (strain, 
polarization, dielectric constant, etc) but also into dynamics and time-correlation, including 
fos instance dielectric loss|4]. Another very interesting possibility of MD is to perform 
thermodynamic integration to obtain free energy differences such as Landau free energies of 
ferroelectricsjl^|. 
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